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ABSTRACT 



A time-dependent, two-dimensional thermodynamic model of an open lead 
in central arctic sea ice is presented. The model is generated by opening 
a lead of finite width and infinite extent in the equilibrium sea-ice cover. 
From this initial condition, the model is integrated numerically as a sea- 
ice cover is reestablished over the lead. The effects of various representa- 
tive advective parameterizations , and temperature and salinity profiles in 
the ocean's surface layer, on the heat flux through the lead and on the 
nature of ice formation are investigated as the lead closes by thermal 
processes alone. Continuity equations involving a horizontal advection 
term and a vertical diffusion term govern heat and salt transport in the 
water. Cooling-induced convective overturn as a mechanism for vertical 
heat and salt transport in the water column is treated through an artifice 
of the vertical diffusion term. 
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I. INTRODUCTION 



A. POLAR AIR-SEA INTERACTIONS 

Heat lost from polar regions together with the heat source of the 
tropics form the basis of the global general circulation. The flux of heat 
into polar regions from lov/er latitudes results from the general meridional 
temperature gradient from polar regions toward the equator due to unequal 
latitudinal heating of the earth's surface. Most of the heat advected into 
polar regions by the atmosphere is lost to space by long wave radiation, 
thus forming a basic mechanism for maintaining planetary radiative 
equilibrium. 

Heat is also advected into polar regions by ocean currents, but due 
to the relative inefficiency of air-sea heat transfer and the relative 
slowness of ocean currents, this transfer is roughly one-tenth that 
attributed to the atmosphere. 

Sea ice is a dominant feature of both arctic and antarctic seas; however 
due to major differences in the physical configuration of sea and land 
masses of these regions, it plays decidedly different roles. The arctic, 
being essentially an ice-covered ocean, is extremely sensitive to climatic 
perturbations since the areal extent and depth of sea ice is critical to the 
regional energy balance. A surplus of about one-third of the annual energy 
influx would remove the arctic ice pack [Maykut and Untersteiner 1969, 1971] 
On the other hand, the antarctic continental mass, covered by 2000 to 4000 
meters of ice, is relatively insensitive to small energy flux perturbations, 
although its sea-ice cover would be affected similarly to that of the arctic 
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Of major import to polar air-sea interactions is the degree to which 
the sea surface is covered by ice. Sea ice acts analogously to a 'lid' 
on the sea, modifying the dynamics and continuity of the sea and the air- 
sea system. Sea-ice cover dramatically reduces the heat transfer between 
ocean and atmosphere by effectively limiting that transfer to a process of 
molecular conduction. Further, by intercepting the vertical flux of 
momentum transferred from the mean wind field to the sea surface by the 
Reynold's stress, it suppresses wind mixing. Finally, it reflects a 
considerable portion of the incoming short-wave radiation, imposes an 
upper limit on the sea's surface-boundary temperature, reduces latent heat 
transfer by impeding evaporation, and, through its changes of phase, serves 
as a buffer, or damping mechanism. Alternatively, heat transfer processes 
are significantly enhanced over open-water areas in proximity to large ice 
masses, particularly during the winter, as surface air temperatures drop 
far below surface water temperatures. Wind and water stresses, in addition 
to directly affecting ocean-atmosphere heat transfer processes, keep the 
ice in nearly continuous motion, thus locally affecting heat transfer as 
well as producing mechanical stresses resulting in the formation of pressure 
ridges, polynya, and open leads. 

Polar air-sea interactions are therefore characterized by complex feed 
back wherein the ocean and atmosphere individually and collectively influence 
the thickness and areal extent of the ice, which in turn has a pronounced 
effect on the exchange in question. 

B. SURVEY OF PREVIOUS RESEARCH 

Range and depth of interest in arctic research has expanded considerably 
since Nansen's epic drift in the FRAM (1893-1896). From pioneering studies 
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until only recently, research emphasis has been primarily of a descriptive 
nature, and it was not until the Soviet establishment of the drifting 
station "North Pole 1 (NP-1)" in 1937 that a systematic study of the 
arctic interior was initiated. Launching of a comparable long-range 
arctic research program by the United States would await the establishment 
of "Fletcher's Ice Island (T-3)" in 1952. 

Attempts to understand the response of sea ice to environmental 
changes has long stimulated the interest of arctic researchers. Quantita- 
tive prediction of these responses has classically involved two primary 
approaches, the empirical and the theoretical. Although providing minimal 
insight to the dynamics of physical processes occurring between air, ice, 
and sea, empirical formulae relating sea-ice growth processes to observed 
surface temperatures have accomplished basic predictive goals. The most 
sophisticated example of this approach is perhaps one advanced by Bilello 
[1961]. Theoretical studies involved analytical solutions and were there- 
fore severely limited in scope and complexity. One of the foremost of 
these analytical treatments is attributed to Kolesnikov [U.S. Naval Civil 
Engineering Laboratory 1966]. 

In the early 1950' s the tenor of arctic research began to reflect 
increased emphasis on the theoretical approach to arctic dynamics as 
comprehensive data from systematic studies became available. It was at 
this point that studies relating the arctic energy balance to the mass 
budget and subsequent research directed toward heat transfer processes 
began to emerge. Development of physical relationships was confined to 
empirical and statistical approaches, or theoretical problems to which 
simple, analytical, approximate solutions could be developed. 
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Classic arctic heat budget studies were presented by several authors. 
Badgley [1961, 1966] summarized and collated available data to deduce a 
typical heat budget at the surface of the central Arctic Ocean. Fletcher 
[1966] discussed annual patterns of atmospheric heat loss for each component 
of the heat budget in relation to the general atmospheric circulation, for 
both ice-covered oceans and ice-free oceans. Vowinckel [Arctic Meteorology 
Research Group 1964] independently investigated the amount of heat conducted 
through the Arctic Ocean ice to check previously calculated values for the 
heat released from the Arctic Ocean. Doronin [1966] presented an analysis 
of heat-balance components of the surface layer in the Arctic Ocean. 

Coachman [1966] proposed a model to explain the observed features of the 
oceanic regime of the Arctic surface layer during the winter. Muench 
[Baffin Bay - North Water Project 1971] discussed the physical oceanography 
of the northern Baffin Bay region in considerable detail. 

It was not until the 1960's and the attendant widened use of the high- 
speed digital computer as a research tool that physical processes could be 
modeled free “of constraints imposed by analytical solutions. Using finite 
difference schemes, a system of integral -differential equations could now 
be solved rapidly within acceptable error limitations. Typical of first- 
generation numerical models in the field is that developed by Untersteiner 
[1966], which predicted thickness and temperature of sea ice, thus allowing 
further insight into the relationship of climatic change to ice production, 
ice decay, and equilibrium ice thickness. The advanced one-dimensional 
thermodynamic model of central Arctic sea ice developed by Maykut and 
Untersteiner [1969, 1971] is a prime example of the present theoretical 
approach. 
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C. PRESENT RESEARCH DIRECTION 



A great deal of research effort is presently being directed toward 
developing numerical models of the atmosphere and of the ocean as well 
as atmosphere-ocean models, with the objective of ultimately producing 
an efficient, reliable global system model . Problems inherent in the 
realization of this objective are numerous and include not only a lack of 
refinement in respective models of the atmosphere and ocean, but also a 
lack of refinement in the technique of coupling the respective models. 
Further, observational input data is in general deficient for these models, 
particularly in the Arctic. Present state-of-the-art examples of models 
of this type are well represented by the sophisticated models of Manabe 
[1969] and Bryan [1969]. 

One of the more ambitious researcfi projects in recent years is the 
Arctic Ice Dynamics Joint Experiment. AIDJEX is essentially a United 
States - Canadian cooperative effort to gain quantitative understanding 
of the interaction between atmospheric sea ice and fluid ocean fields of 
motion. In progress at the time of this writing, AIDJEX is unique in its 
observational approach since, in order to gain a full understanding of 
interacting fields, it involves the collaboration of many research groups 
in making coordinated measurements of interacting fields of motion, stress, 
and strain - observing appropriate time and space scales - over a minimum 
period of one year. The fruits of this ambitious and imaginative study, 
in terms of fresh insights, controlled data, and international cooperation 
appear to be vast already. 
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D. CONCEPTUALLY-RELATED STUDIES 



Studies of Arctic dynamics, energy balances and mass budgets range 
over a wide spectrum. Scientific literature reflects numerous experimental 
and theoretical models of microstructure and micro-scale energy transfer 
processes between ocean and atmosphere, ice and atmosphere, and ice and 
ocean, essentially occurring in boundary-layer regions. Pertinent to the 
work presented herein is research performed by Foster [1968, 1969] in which 
haline convection induced by the freezing of sea water was studied, and by 
Lake and Lewis [1970] who studied salt rejection by sea ice during growth. 
Work performed by Lewis and Walker [1970] illustrates seasonal changes in 
the temperature and salinity profiles beneath annual sea ice cover, and 
Stern and Turner [1969] reported on the formation of salt fingers and 
convecting layers. 



Directly related to the model proposed herein is the macroscale advanced 
model derived by Maykut and Untersteiner [1969, 1971] to forecast ice 
temperature and thickness, which defines the role of each component of the 
energy budget in relation to its effects upon the ice. Although comprehen- 
sive in its treatment of physical processes related to components of the 
energy budget, the model admittedly has several shortcomings. Physical 
limitations are due to uncertainties in environmental data and include 
the necessary assumption of a constant oceanic heat flux and generally 
course albedo and energy-flux data. Limitations, arising from an inability 
to model certain processes include: 1) the neglect of mechanical stresses 

produced by winds and surface currents, 2) treatment of turbulent fluxes 
at boundaries as being independent of growth rates or the physical states 
of ice and snow, 3) imprecise treatment of downward heat transport from 
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melting snow, 4) undefined shapes of salinity profiles for various ice 
conditions, and 5) neglect of heat storage by melt ponds. Leads were 
considered only so far as they indirectly affected the energy fluxes over 
the bare ice. In spite of its limitations, this will certainly be a 
cornerstone for future thermodynamic ice models and provides a number of 
conceptual bridges used in the present work. 

E. THESIS OBJECTIVE 

It . is known that significant heat transfer occurs between the atmosphere 
and the ocean through open leads in the Arctic sea ice cover. As found 
experimentally by Badgley [1961], sensible heat loss is at least two orders 
of magnitude greater from open leads than from perennial sea ice. 

With the need for further refinement in models of the ocean and the 
atmosphere, as well as in model coupling schem.es, the requirement for 
mesoscale examination of ocean-ice-atmosphere heat- transfer processes is 
immediate and apparent. With appropriate observational input data becoming 
increasingly- more abundant, informative numerical models of mesoscale 
processes appear promising for furthering the understanding of physical 
processes and for realizing control of mass-energy fluxes in the arctic. 

The objective of this study was to develop a thermodynamic model of an 
open lead in central arctic sea-ice in v;hich hydrodynamic processes of 
convection, advection, and diffusion of heat and salt were represented. 

By integrating the model numerically, the effects of various representative 
advective parametizations , and temperature and salinity profiles in the 
ocean's surface layer, on the heat flux through the lead and on the nature 
of ice formation, could be investigated as the lead reclosed by thermal 
processes alone. 
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II. NATURE OF THE PROBLEM 



Since sea ice essentially insulates the ocean from the atmosphere and 
suppresses important thermodynamic processes, it, along with heat-advecting 
ocean currents, becomes a primary factor in determining the vertical 
structure in Arctic seas. However, the dynamic transition from ice-free to 
ice-covered seas is not fully understood. In particular, the nature of 
the fluxes of salt and heat, as well as the structure and depth of unstable 
convective overturn during initial phases of the freezing process, require 
study. 



A. DYNAMICS OF THE SEA WATER FREEZING PROCESS 

The structure of the active surface layer of the Arctic Ocean is 
generally .horizontally uniform with a very stable density gradient. How- 
ever, during the freezing process the density gradient may be significantly 
altered by three primary mechanism: 1) progressive cooling of surface 

waters from above, 2) brine expulsion by the ice during the freezing pro- 
cess, and 3) advection of heat and/or salt by surface currents. 

The fundamental physical relation which makes the sea water freezing 
process unique may be understood by examining the curves of maximum density 
and freezing point for sea v/ater, plotted with respect to temperature versus 
chlorinity or salinity (Figure 1). Specifically, for waters of salinity 
greater than 24.70 parts per thousand, the temperature of maximum density of 
the liquid state will always be lower than the liquid's freezing-point 
temperature under pressures existing in nature. Any process that tends to 
produce a vertical density instability will induce convective overturn in 
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Figure 1. 

Temperatures of Sea Water 
Maximum Density and 
Freezing-Point 




Chlorinity o/oo 
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the water column, resulting in mixing of the unstable segment such as 
would be produced in the freezing process by surface cooling and freezing- 
associated brine expulsion. Convective overturn will therefore be continuous 
in the unstable segment of the liquid water column from the onset of surface 
cooling to the point of total change-of-phase of the water column. A some- 
what less obvious but straightforward point should be mentioned in qualifi- 
cation of the foregoing. From an arbitrary water temperature above the 
freezing point, down to, but not including the freezing point, convection 
is temperature controlled during cooling since there is no ice formed, 
hence, no brine expulsion. During cooling of surface waters at the freezing 
point for existing surface salinities, convection is primarily salinity 
controlled since brine is being continuously expelled from newly forming 
ice. Within this latter water temperature regime, temperature changes play 
a minor role in effecting changes in surface water densities. 

As convective overturn is a mechanism for inhibiting freezing by 
continuously transferring heat from depth to the active boundary at which 
freezing is “taking place, so it is also a mechanism for effectively removing 
heat from depths significantly greater than could be effected by molecular 
diffusion. An equilibrium ice thickness is reached when the resultant heat 
flux upv;ard through convection, advection and diffusion in the water column 
balances the upward resultant heat flux through the ice, at the ice-water 
boundary. Prior to this point, ice is accreted in depth since the upward 
heat flux through the ice exceeds the upward heat flux in the v/ater column 
arriving at the ice-water boundary. 
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B. DISTRIBUTION OF VARIABLES 

The distribution of temperature, salinity and velocity are considered 
within a volume of water under a refreezing lead. The water is assumed to 
be inviscid, horizontally homogeneous, horizontally isotropic, irrotational , 
and incompressible. Accelerations resulting from the coriolis effect and 
horizontal pressure gradients are assumed to be negligible. Momentum, heat 
energy, and salt mass are conserved within the volume such that no velocity, 
temperature, or salinity sources or sinks exist except at the volume's 
boundaries. Furthermore, currents are assumed to be steady-state and 
uniform both horizontally and vertically. 

Using the classic box-model approach, these conservative fields are 
represented continuously within the volume by individual continuity 
equations. For the vector field of velocity, where K is the diffusivity 
and S,, is a source term, 

V 

= V-KW+ (1) 

which reduces to 

v(t,x,y,z) = V = constant (2) 

i.e., the velocity is specified as uniform and steady. 

Similarly, the scalar fields of heat and salt, here collectively 
depicted by the scalar quantity 0, are represented by 

D0/ot = V . K V 0 + (3) 

where K is the diffusivity and is a source term. 
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This reduces to 



■ V - v0 + . K(z) v^0 (4) 

where is the horizontal diffusivity and is the vertical del operator 
since , 



K(x) = K(y) = = constant 

The local, or time-rate-of-change , is therefore equal to the sum of the field- 
rate-of-change and the diffusive-rate-of-change. 



C. HYDRODYNAMIC PROCESSES 

Mechanisms for the hydrodynamic transport of heat and salt within the 
lead are advection and diffusion. However, from the dynamics of the sea- 
water freezing process, it is knov/n that convection plays a major role in 
transporting both heat and salt in the water column. Although it is not 
represented by a discrete term in the foregoing continuity equation, the 
convective transport term may be thought of as simply a temporary (when 
vertical density instabilities arise) modification of the vertical diffusion 
term taking the form 



K(z)-oo 

which in essence results in free vertical (convective) transport. 



(5) 



D. MASS AND ENERGY BALANCE 

With the foregoing analyses of continuity relations and mechanisms of 
hydrodynamic transport existing within the box model of the open lead, the 
mass and energy balances existing at the volume's boundaries may be 
examined. 
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Over an arbitrary surface area of a reclosing lead, a given change of 
heat, aQ, at the water surface must equal the heat added by advection and 
diffusion at upstream boundaries, Q^, plus the heat added from depth in the 
lead by convection, advection, and diffusion, Q^, plus the heat added by 
surface ice accretion, Q^, plus the open-water heat loss at the surface by 
conduction, evaporation and radiation, plus the heat loss through 

the ice cover by conduction, , minus the heat diffused across vertical 
boundaries normal to advective transport, Q^, minus the heat transported 
across the downstream boundaries by advection and diffusion, Q^j, such that 
aQ = Qu Qz Qf ^{Qg+Qg+Q^) + Q^(l-W) - Qp - Qd (6) 

Where W is the fraction of open water comprising the surface of the lead at 
any given time. The heat added from depth in the lead is heat ultimately 
derived from the Bering Sea Water (in the Canadian Basin) and from the 
Atlantic layer. 

By similar analysis, a given change of salinity, aS, at the water's 
surface must equal the salt added by advection and diffusion at upstream 
boundaries, Sjj, plus the salt added from depth in the lead by convection, 
advection, and diffusion, S^, plus the salt added by surface ice accretion 
(salt rejection), S^, plus the salt lost through surface ice melting (salt 
dilution), plus the salt added by surface evaporation, Sg, minus the 
salt diffused across vertical boundaries normal to advective transport, 

S^, minus the salt transported across downstream boundaries by advection and 
diffusion, S^, such that 

AS = + (l-W)(Sf+S^) + W(Sg) - (7) 

where again W is the fraction of open water comprising the surface of the 
lead at any given time. The salt added from depth in the lead is derived 
from the generally strong positive halocline existing at depths below 25 to 
50 meters. 
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Qualitatively then, before the equilibrium sea ice cover has been 
reestablished, there will exist differentials in the temperature and 
salinity profiles of water columns. A primary process of interest in the 
refreezing lead is that of surface ice formation, which lends a unique 
character to lead hydrodynamics. The heat added/lost by surface ice 
accretion/melting is a function of the incremental change in the mass of 
ice,AI, and the heat of formation of sea ice, , such that 

Qf = -Qm = Li*AI (8) 

Coupled directly to this process is the process of brine expulsion (salt 
rejection) by the sea ice during growth. The salt added/lost by surface ice 
accretion/melting is a function of the incremental change in the mass of ice, 
Al, the salinity of the sea ice, , (which in turn is dependent upon rate- 
of-freezing, ice temperature, and age of the ice), and the salinity of the 
parent sea water, S.^y, such that 

Sf = -Sm = f( Al,S^-,S,y) (9) 

Qualitative analysis of the salt rejection process by constituent salts is 
not appropos to either the scale or the nature of the study at hand; there- 
fore salt rejection was treated simply as a process of water-salinity 
modification. 
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III. EXPERIMENTAL PROCEDURE 



A. FORMULATION OF THE MODEL 

1 . Establishment of a Representative Cross Section 

In deciding upon what physical dimensions would be appropriate for 
the model, it appeared most suitable to represent a small open lead in 
central arctic sea ice such as one would expect to observe between early 
to late winter. Therefore, the width of the lead was chosen to be 150 
meters, with an additional 50 meters under equilibrium sea-ice cover on 
the downstream side of the model lead to facilitate studying resultant 
thermodynamics and hydrodynami cs in that region. 

The model's depth extent was set at 50 meters based upon lower 
boundary condition considerations. Over much of the Arctic Basin the 
vertical temperature gradient at that depth is zero or very slightly posi- 
tive. Since a no-flux lower boundary condition was desirable in order to 
better study the effects of advection of conservative constituents in the 
lead, that depth was chosen. 

It was felt that the model could be simplified without the loss of 
any essential physics if derivatives along the lead (Y-axis) v;ere set to 
zero. This reduced the problem to one of looking at a cross-sectional 
distribution in the lead. Figure 2 is a pictorial depiction of the open 
lead. 

2. Adaptation of the Continuity Equations to the Model 

Using the simplification of the previous paragraph which allows the 
use of a cross-sectional area rather than a true volume to represent the 
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Figure 2. 

Pictorial Depiction of Open Lead 
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lead, the horizontal advection and horizontal diffusion terms describing 
changes along the Y-axis were eliminated from equation (4). 

The problem was further simplified by assuming that all horizontal 
diffusion was negligible relative to horizontal advection, and that verti- 
cal advection could be disregarded. 

Applying the foregoing simplifications to the expanded form of 
equation (4), and choosing the state variables of temperature and salinity 
as the relevant dependent variables, the final form of the equations used 
to describe the conservative scalar fields of temperature and salinity 
within the lead was 



0=0 (t,x,z) 

■ 3. Initial Conditions 

Sea-water temperatures and salinities within the lead at the 
instant that'the lead is generated are represented by local stable tem- 
perature and salinity profiles previously existing beneath the equilibrium 
sea-ice cover. The temperature and salinity fields are both horizontally 
uniform and vary only with depth in the lead such that 



The surface of the newly-generated lead is initially taken to be 
ice-free over the left (upstream) 150 meters, and ice-covered with equili- 
grium thickness ice over the right (down-stream) 50 meters. Surface 
temperatures are initially at the freezing point, in temperature equilibrium 
with the preexisting sea-ice cover. 




( 10 ) 



where , 



T(0,X,Z) = T(Z) 
S(0,X,Z) = S(Z) 



( 11 ) 

(12) 
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4. Boundary Conditions 



The model's left, or upstream boundary, may be regarded as a 
source boundary across which sea water having a given temperature and 
salinity profile flows into the lead. Sea-water temperature and salinity 
profiles are taken as steady-state functions of depth alone, identical to 
the lead's initial temperature and salinity profiles such that 



energy fluxes and is therefore the most dynamically active boundary. The 
boundary condition on temperature across the entire surface is that of an 
upward (or downward) diffusive flux of heat across the water surface such 
that 



where aT is the air-water temperature difference (an externally specified 
parameter), is the surface ice thickness, S^- is the surface ice salinity 
(an externally specified constant parameter), and and K-j are the co- 
efficients of thermal eddy diffusivity of sea v/ater and thermal conducti- 
vity of sea ice respectively. When no ice is present at the surface, the 
surface heat flux is taken as the cumulative open-water sensible, latent, 
and radiative heat fluxes, which are specified in the model as external 
parameters according to the air-water temperature difference. Since the 
model's surface at all times represents the water surface, the fluxes at 
that boundary are matched as ice is formed, such that the flux through the 
water surface is always equal to the flux through the ice to the atmosphere, 
but at no time does it exceed the maximum open-water flux to the atmosphere. 



T(t,0,Z) = T(Z) 
S(t,0,Z) = S(Z) 



(13) 

(14) 



The model's surface boundary is the primary region of mass and 




(15) 



25 



The boundary condition on sea water salinity at the surface 
specifies a no-diffusive flux condition betv/een air and water as well as 
between ice and water and ice and air such that 

S/g_ S(t,x,o) = 0 (16) 

However, salt rejection from ice to water during the freezing process (as 
well as salt dilution during melting) is represented as a discrete process 
of physically adding the change in salinity produced by sea-ice accretion 
or melting to the surface waters such that 

S(t,x,o)=» S(t,x,o) + f ( Al , aS) (17) 

where AI is the incremental change in the mass of sea ice, and AS is 
the salinity differential existing between the sea ice and its parent sea 
water. 

The bottom boundary, as previously discussed, is characterized 



as a no-flux boundary, thus 






T(t,X,h) 


= 0 


(18) 


S(t,X,h) 


= 0 


(19) 



where h is the model's depth extent. 

The right, or downstream boundary, is also characterized as a 



no-diffusive flux boundary where whatever temperature-salinity field is 
generated by the thermodynamic processes of the lead is simply advected 
out. However, as mentioned previously, the right 50 meters of the lead 
model's surface is overlain by an equilibrium sea-ice cover. This boundary 
condition is therefore representative of a return towards the equilibrium 
vertical profiles. At the model's right boundary then 
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ax T(t,W,Z) = 0 



( 20 ) 



3/ 



S(t,W,Z) =0 (21) 

where W is the model's width. 

Figure 3 describes processes considered in the model cross section. 
Figure 4 schematically represents the model cross section with its governing 
continuity equations and boundary conditions. In this figure, the surface 
boundary condition is simply a statement of the balance of fluxes at the 
ice-v/ater boundary. With no ice cover (5=0), the upward heat flux through 
the v/ater surface, F^. is equal to the sum of open-water sensible (F^), latent 
(Fg), and radiative (F^) heat fluxes. When an ice cover is present (5>0), the 
upward heat flux through the water surface (z =0‘'') is equal to the upv/ard heat 
flux through the ice (^ =0") minus the heat added by surface ice accretion. 

The above relations are presented in terms of temperature gradients in Figure 
4 as they are applied to the governing continuity equation for temperature. 
Symbols not previously defined are volume density of sea ice; q^ , heat of 

formation of sea ice; 5 ice thickness (volume of ice per unit surface area). 

5. The Finite-Difference Scheme 

A purely explicit finite difference scheme was selected to numeri- 
cally represent the linear parabolic partial differential equation (10). 

The time rate-of-change term was represented by forward differencing, the 
field rate-of-change term by backward (upstream) differencing, and the 
diffusive rate-of change term by central differencing. Thus over the interior 
(non-boundary) region of the model 



^ t 



A X 



(AZ)2 



( 22 ) 
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Figure 3. 

Description of Processes Considered in Model 



( 2 ) 




(1) Advection of steady-state temperature-salinity field into the 
lead from under equiliorium ice pack. 

(2) Heat loss to the atmosphere as a function of ice thickness, 
ice accretion/melting, and salt rejection/dilution. 

(3) Heat and salt transport by convective overturn, advection and 
eddy diffusion. 

(4) No transport across boundary. 

(5) Advection of temperature-salinity field generated by lead 
thermodynamics out of the lead model. 
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Figure 4. . . 

Model Cross Section with Governing Continuity Equations and Boundary Conditions. 
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where , 




o"''* t 

JAX, kAZ 



(23) 



At the left boundary, equation (22) was inapplicable since tempera- 
ture and salinity v/ere both steady-state functions of depth alone. There- 
fore, at this boundary at all times 



0 



n+1 

j,k 




(24) 



At the surface boundary, the diffusive term for temperature was 
modified in equation (22) since 



or. 



‘j,k+l 'j,k-l 
2 A Z 



— A T / r7 T 



= ^ T 

V ' 



(25) 



(26) 



The diffusive salt term at the surface was modified in equation (22) 



since 



^j,k+l 



^j,k-l 



2 A I 



(27) 



or. 



c-n _ .n 
^j,k-l ^j,k+l 



(28) 



At the bottom boundary, the diffusive rate-of-change terms of both 
temperature and salinity were modified in equation (22) since 



2AZ 



= 0 



(29) 
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or, 




( 30 ) 



Finally, at the right boundary, equation (22) remains unmodified 
since advection is the only transport mechanism operating across it. 

, Stability analysis of the finite difference equation using the 
Fourier series method [Smith 1965] shows that 



where < represents an arbitrary error introduced into the problem. There- 
fore, the finite-difference equation (22) is stable when the absolute value 
of the amplification factor remains less than or equal to unity 



6. Representation of Physical Processes 

Convective overturn is treated nominally through an artifice of the 
vertical diffusion term as in equation (5). Numerically however, convective 
overturn and its associated mixing is accomplished through arithmetic 
averaging of temperatures and salinities within regions of vertical density 
instability as determined by sigma-t values, followed by recomputation of 
the resultant sigma-t value for the newly mixed (and neutrally stable) 
region. Subroutine OVRTRN accomplishes convective overturn in the model. 

Advection of heat and salt (temperature and salinity) in the lead is 
driven by a current, u, in the positive-x direction. Values and depth pro- 
files of the current as used in the model are given in subparagraph B of this 
section. 




(31) 




(32) 
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Diffusion of heat and salt in the vertical are determined by the 
coefficients of vertical temperature and salinity diffusivities , both set 

p _i 

at constant values of 10 cm sec in the model. 

Formation of sea ice is accomplished as in equation (8). More 
specifically however, when a given mass of water, Aw, is cooled AT degrees 
centigrade below its salinity dependent freezing-point temperature, the 
incremental mass of ice formed, A I, is 

A I = S(^W.AT)/, (33) 

where Cp is the specific heat of sea water and L^- is the heat of formation 
of sea ice. 

Salt rejection (or dilution) resulting from surface sea ice accretion 
(or melting) is described in equation (17). Specifically, if A I is the 
incremental change in the mass of sea ice and AS is the salinity differential 
existing between a unit mass of parent sea water and a unit mass of its 
generated sea ice, then the incremental salinity change, ASp, in the surface 
waters resulting from that freezing or melting process is 

A Sf^ = A I-AS (34) 

7. Physical Quantities and Relations 

Sea water chlorinity is related to salinity [Knudsen 1902] by the 

relation 

C1=(S - 0.030)71.8050 (35) 

Freezing-point depression is empirically related to sea water 

chlorinity by Thompson [Sverdrup, Johnson and Fleming 1942] by the relation 

T^ = -0.0966-Cl - 0.0000052 -Ci3 (36) 

fP 
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The density of sea water as a function of salinity and temperature 



(sigma-t) [Knudsen 1902] is empirically computed through 



<r = 
0 


»l(Cl) 


(37) 




^2(T) 


(38) 


®T = 


^3(T) 


(39) 


^T = 


^4(T) 


(40) 




^5^^T’ *^o’ ^T’ ®T^ 


(41) 



The specific heat at constant pressure of sea water is empirically 
related to its salinity at normal temperature and pressure by Kuwahara 
[Sverdrup, Johnson and Fleming 1942] by the relation 

Cp = 1.005 - 0.004136-S + 0.0001098‘S2 

- 0.000001 324- $3 (42) 

The thermal conductivity of sea ice was determined using the relation- 
ship established by Unterstcinsr [Maykut and Untersteiner !959, 1971]. 

k. = , T<0 (43) 

where, 

the subscript f refers to pure ice, 
k^ = 0.00486 cal/cm-sec-°C, 

p 

b = 0.28 cal-cm /g-sec, 

T = temperature of ice in degrees Celsius, 

O 

S(z) = ice salinity at depth z in g/cm . 

The heat of formation of sea ice was represented by the simplified 
relationship [Pounder 1965] 

L. = L^.(1--|-) (44) 

where, 

= 79.77 cal/g, 

= salinity of sea ice, g/kg, 

S = salinity of parent sea water, g/kg. 
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B. ENVIRONMENTAL DATA INPUTS 

Environmental data used in the model fell into two basic categories: 
that which remained constant throughout all runs of the model, and that 
which v/as variable from run-to-run. Representative of the former were 
initial equilibrium-ice thickness, set at 200 cm, an ice salinity of 8.0 
g/kg, and an ice density of 0.91 g/cm^. Further, representative open water 
sensible, latent, and radiative heat flux densities as presented by Muench 
[Baffin Bay-North Water Project 1971] were applied in the model. However, 
as ice formed at the surface, the latent heat flux density was modified so 
as to decrease exponentially with increasing ice thickness. The radiative 
heat flux density v;as modified so as to decrease from open water values as 
ice formed, approaching a constant value representative of the radiative 
flux density over pack ice as also presented by Muench. The sensible heat 
flux density v/as modified internally by the relation shown in Figure 4 as 
surface ice thicknesses increased. All surface heat flux densities were 
representative of an air-water temperature difference of 25°C with no wind. 
Finally, the coefficients of vertical temperature and salinity eddy diffusi- 

p 

vity were set at 10.0 cm / sec. 

Primary variable parameters in the model were horizontal current magni- 
tudes, vertical temperature profiles, and vertical salinity profiles. 
Currents applied in the model were within the range of zero to 10.0 cm/sec. 
Extreme test case values were 0.0 and 10.0 cm/sec. Additional test case 
values of 2.0 and 4.0 cm/sec were used to study one particular temperature- 
salinity field combination with its associated hydrodynamics in further 
detail. A value of 7.0 cm/sec was chosen to represent a typical horizontal 
current occurring in the lead, and was subsequently used as the real case 
value. 
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Test case values of vertical temperature and salinity profiles 
representative of extremes in vertical gradient were applied. Profiles 
A represented strong gradients, while profiles B represented weak 
gradients (Table I). Real case temperature and salinity profiles, 
qualitatively representative of early winter (profile C) and late 
winter (profile D) conditions in the central arctic [Coachman and 
Barnes 1961, 1962], were chosen to represent the typical lead. 

Table I summarizes the profile designations of the three variable 
parameters and their various values. 



C. MODEL INTEGRATION PLAN 

In order to gain maximum insight to physical processes occurring in 
the lead, as well as to study the effects of current magnitude, vertical 
temperature gradient, and vertical salinity gradient on lead thermodynamics, 
the series of test variables vias run in tns model, issprescnted by runt. 



one through eight, these test cases involved running all possible combina- 
tions of current magnitude (u) profiles A and B, temperature (T) profiles 
A and B, and salinity (S) profiles A and B. These test cases were 
integrated over a lead reclosure period of tv70 days. 

As a result of the interesting dynamics occurring in run 6, additional 
runs 9, 10, and 14 were scheduled using temperature and salinity profiles 
identical to run 6 (and 2), but with current magnitude profiles C, D, and 
F respectively, in order to further investigate these phenomena. Integra- 
tion was again for a period of two days. 

Finally, the real-case variables, representative of naturally 
occurring parameters, were put into the model and integrated over two-day 
periods in runs 11 and 12, and over a 30 day period in run 13. 



I 
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Table I. Environmental Variable Profiles 



Horizontal Current Magnitudes (cm/sec) 



Depth 

(m) 


Profile 


A 


B 


C 


D 


E 


F 


0-50 


10.0 


1.0 


2.0 


4.0 


7.0 


0.0 



Water Temperatures (OC) 



Depth 

(m) 


Profile 


A 


B 


C 


• D 


0 


O 

LO 

1 


-1.50 


-1.68 


-1 .82 


5 


-1 .40 


-1.49 


-1.68 


-1.82 


10 


-1.30 


-1.48 


-1.68 


00 

ro 


15 


-1 .20 


-1.47 


-1.67 


-1.82 


20 


-1.10 


r— 

1 


-1.66 


-1.81 


25 


-1.00 


-1 .45 


1 

cn 


-1.80 


30 


-0.90 


-1.44 


-1 .64 


-1.79 


35 


-0.80 


-1.43 


-1.63 


1 

00 


40 


1 

O 

o 


-1.42 


-1.62 


1 " "7 *7 

1 


45 


-0.60 


-1 .41 


-1 .61 


-1.76 


50 


-0.50 


-1.40 


-1.60 


-1.75 



Water Salinities (g/kg) 



Depth 

(m) 


Profi le 


A 


B 


C 


D 


0 


28.00 


28.00 


31.00 


33.51 


5 


28.50 


28.01 


31.01 


33.52 


10 


29.00 


28.02 


31.02 


33.53 


15 


29.50 


28.03 


31.06 


33.55 


20 


30.00 


28.04 


31.15 


33.57 


25 


30.50 


28.05 


31.25 


33.60 


30 


31 .00 


28.06 


31 .35 


33.63 


35 


31 .50 


28.07 


31 .45 


33.66 


40 


32.00 


28.08 


31.55 


33.69 


45 


32.50 


28.09 


31 .65 


33.72 


50 


33.00 


28.10 


31.75 


33.75 
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Table II presents a summary of runs conducted with their respective 
variable combinations. 



D. MODEL LIMITATIONS 



In addition to the assumptions made earlier in arriving at the 
simplified model of an open lead, basic limitations in the model s 
capabilities are: 1) the atmosphere is treated as an external parameter 

such that it affects lead dynamics but cannot be affected in turn by the 
lead; 2) sea ice is treated as a quasiexternal parameter since it is 
depicted as a simple, uniform mass floating on the water's surface. 

Density and salinity are constant throughout this mass in both space and time. 



and the vertical temperature gradient is taken to be linear throughout. Its 
thermal conductivity is determined only at the ice-water boundary, which is 
warmest and therefore limiting through the relation given in equation (43), 
and an upper limit is applied to the upward heat-flux density through the 
ice such that at no time is the flux through the ice allowed to exceed that 



over open water; 3) mechanical lead closure (such as by wind-drifted ice) 
is disallowed; 4) snow cover over the lead's surface is not modelled; 5) 
super-cooling of surface waters [Coachman 1966] is not portrayed; 6) as 
convection-induced mixing occurs instantaneously in the model, so also 
does surface ice melting whenever heat is brought to the surface by 
convection. This results in an artificiality in that, although heat and 



salt are conserved in the model, heat brought to the surface under an ice 
cover immediately melts ice until no further surface heat is available for 
melting. Therefore, under this limitation, it is not possible for surface 
heat to be advected downstream under an ice cover to promote downstream 
melting. Only in cases when no ice is present to be melted may this heat 
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Table II. Schedule of Runs and Variable Combinations 



Run 


Variable Profile 


Period of Integration 
(days) 


U 


T 


s 


1 


A 


A 


A 


2 


2 


A 


A 


B 


2 


3 


A 


B 


A 


2 


4 


A 


B 


B 


2 


5 


B 


A 


A 


2 


6 


B 


A 


B 


2 


7 


B 


B 


A 


2 


8 


B 


B 


B 


2 


9 


C 


A 


B 


2 


10 


D 


A 


B 


2 


11 


E 


C 


C 


2 


12 


E 


D 


D 


2 


13 


E 


n 


n 


30 


14 


F 


M 


B 


2 
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be advected downstream; 7) in situations wherein vertical salinity gradients 
existing in the water column are weaker than a certain minimum gradient, 
the model will over-mix the water column in depth. Specifically, the model 
forms ice in response to surface cooling, thereby rejecting salt into the 
surface water layer of depth Z. The model then mixes the induced un- 
stable region as determined by relative sigma-t values of density. If 
however, there is heat available at some shallower depth than that to 
which the water column was mixed which could have decreased the initial 
amount of ice formed, then although salt, ice, and heat are ultimately 
conserved, potential energy is not conserved. Initial freezing and resul- 
tant salt rejection have induced convection to a greater depth than would 
have actually been the case in a continuous process. Therefore, for 
vertical salinity gradients weaker than a certain minimum gradient, the 
model's numerical convective process acts as a potential energy source. 
Quantitatively, the limiting minimum vertical salinity gradient which can 
be accurately treated by this model, in g/kg per vertical grid space, is 

. (S^, - sp • ai/A Z (45) 

where is the ice density, is the salinity of sea water, S^. is the 

salinity of sea ice, A I is the incremental change in the volume of ice, 
and AZ is the vertical grid spacing. In this model, the limiting gradient 
is in the neighborhood of 1 X 10"^ g/kg per vertical grid space, v/hich is 
lower by an order of magnitude than the weakest salinity gradient used. 
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IV. PRESENTATION OF DATA 



Some results of integrating the model using the combinations of 
environmental variables summarized in sections III.B. and III.C. are 
presented in tabular and graphic form. Three primary output variables are 
treated: surface ice thickness, surface heat loss to the atmosphere, and 

depth of penetration of surface-freezing induced convective overturn. 
Throughout the section, these three output variables are related by run 
number (section III.C.) to the input variables of current magnitude, 
vertical temperature profile and vertical salinity profile (section III.B.). 

A. TEST CASES 

Ice thicknesses are treated in two ways: by time growth rate and by 

spatial profile across the lead. Values for growth rates were arrived at 
by tabulating the maximum ice thickness occurring anywhere in the lead as 
a function of. time. Values for spatial profiles were taken as the ice 
thicknesses across the model's surface (200 meters) after 48 hours of 
freezing. Values over the right 50 meters of the model's surface, initially 
covered by 200 cm of ice, are tabulated to show the downstream surface 
effects of lead thermodynamics on the equilibrium ice cover. 

Surface heat losses to the atmosphere are presented as cumulative heat- 
loss densities (heat loss per unit surface area) occurring at each grid 
point over the surface of the initially open lead after 48 hours of lead 
reclosure. Additionally, these surface heat-loss densities are averaged and 
tabulated separately as mean surface heat-loss densities. 
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Finally, maximum depths of penetration of surface-freezing induced 
convective overturn are tabulated. Due to the numerical nature of the 
model however, these depths should be treated as bounding depths rather 
than precise depths as would be indicated for continuous schemes. 

B. REAL CASES 

Ice thicknesses, surface heat losses to the atmosphere, and depths of 
penetration of surface- freezing induced convective overturn are treated in 
an analogous manner. Additionally however, ice thickness and mean surface 
heat-loss densities are presented as functions of time for run number 13, 
in which the model was integrated to represent 30 days of lead reclosure. 
Values for surface ice thicknesses were determined analogously to those of 
48-hour ice-growth rates. 
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Ice Growth Rates, Test Cases. 
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Figurs 5. Ic6 Growth Ratos, Test-Caso Runs 
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Time After Lead Generation (hrs) 
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Time After Lead Generation (hrs) 



Thicknesses (cm) Across Lead Surface 4B Hours After Lead Generation 
— Run Number r— 
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Figure 7. Ice Thickness Profiles, Test-Case Runs 
* from Table IV 
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Distance Downstream of Lead's Source Boundary (m) 



Figure 8. Ice Thickness Profiles, Test-Case Runs 
* from Table IV 
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Distance Downstream of Lead's Source Boundary (m) 
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Table V. Surface Heat-Loss Densities, Test Cases 
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Figure 9. Surface Heat-Loss Densities, Test-Case Runs 
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Distance Downstream of Lead's Source Boundary (m) 



Figure 10. Surface Heat-Loss Densities, Test-Case Runs 
* from Table V 
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Table VI. Mean Surface Heat-Loss Densities, Test Cases. 



48-Hour Mean-Surface Heat Losses 
(x 102 ca1/cm^) 



Run 

Number 


AQ 


1 


14.9 


2 


17.0 


3 


13.7 


4 


14.9 


5 


14.9 


6 


23.2 


7 


13.7 


8 


14.5 



Table VII. Maximum Depth of Penetration of Convective Overturn, 
Test Cases. 



Run 

Number 


Deptli 

(m) 


1 


5 


2 


10 


3 


5 


4 


10 


5 


5 


6 


20 


7 


5 


8 


20 
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Table VIII. Ice Growth Rates, Real Cases. 
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Figure 11. Ice Growth Rates, Real Cases. 
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Table IX. Ice Thickness Profiles, Real Cases. 



Ice Thicknesses (cm) Across Lead Surface 
48 Hours After Lead Generation. 



X 


Run Number 


(m) 


11 


12 


10 


24.7 


24.6 


20 


24.7 


24.6 


30 


24.7 


24.6 


40 


24.7 


24.6 


50 


24.6 


24.3 


60 


24.3 


24.0 


70 


24.4 


24.2 


80 


24.4 


24.2 


90 


24.4 


24.1 


100 


24.4 


23.9 


no 


24.3 


23.9 


120 


24.2 

1 


! 23.6 i 

1 


130 


24.0 


23.5 


140 


24.0 


23.4 


150 


203. 


203. 


160 


203. 


203. 


170 


203.- 


203. 


180 


203. 


203. 


190 


203. 


203. 


200 


203. 


203. 
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Figure 12. Ice Thickness Profiles, Real Cases 
* from Table IX 
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Distance Downstream of Lead's Source Boundary (m) 
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Table X. Surface Heat-Loss Densities, Real Cases. 



48-Hour Surface Heat Losses 
(x 10^ cal/cm^) 



X 


Run Number 


(m) 


11 


12 


10 


13.6 


13.9 


20 


13.6 


13.9 


30 


13.6 


13.9 


40 


13.6 


13.9 


50 


13.7 


14.3 


60 


14.2 


14.7 


70 


14.0 


14.4 


80 


13.9 


14.3 


90 


13.8 


14.3 


100 


13.8 


14.3 


110 


13.8 


14.3 


120 


13.8 


14.3 


130 


13.8 


14.3 


140 


13.8 


14.2 


150 


1.99 


2.01 


160 


1.99 


2.01 


170 


1.99 


2.01 


180 


1.99 


2.01 


190 


1 .99 


2.01 


200 


1.99 


2.01 






Figure 13. Surface Heat-Loss Densities, Real Cases. 
* from Table X 
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Distance Dov/ristream of Lead's Source Boundary (m. ) 
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Table XI. Mean Surface Heat-Loss Densities, Real Cases. 



48-Hour Mean Surface Heat Losses (x 10^ cal/cm^) 



Run 

Number 


aQ 


11 


13.8 


12 


14.2 



Table XII. Maximum Depth of Penetration of Convective Overturn, 
Real Cases. 



Run 


Depth 


Number 


(m. ) 


11 


15 


12 


15 
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Table XIII. Ice Growth Rate, Run 13. 



Maximum Ice Thickness (cm) Occurring in the 
Lead's Surface at Five-Day Intervals. 



Day 


Ice 

■ Thickness 


5 


40.6 


10 


59.6 


15 


74.9 


20 


88.3 


25 


101. 


30 


112. 



Table XIV. Mean Surface Heat-Loss Densities, Run 13. 

Mean Surface Heat-Loss Densities (xlO^ cal/cm^) 
at Five-Day Intervals. 



Day 


aQ 


5 


23.5 


10 


34.7 


15 


43.7 


20 


51.6 


25 


58.8 


30 


65.5 
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Figure 14. Ice Growth Rate, Run 13. 
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Figure 15. Mean Surface Heat-Loss Densities, Run 13. 
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V. CONCLUSIONS 



Through the preceding data and the time-dependent temperature-salinity- 
density fields generated in the lead model, numerous insights to the process 
of lead reclosure v/ere afforded, as v/ell as to the hydrodynamic validity of 
the model under certain conditions. 

Perhaps most pronounced v/as the very strong role played by the near- 
surface vertical salinity gradient in maintaining water-column stability. 

A strong halocline near the surface tended to insulate deeper waters in the 
lead from the thermodynamic activity occurring at the surface. Surface 
waters whose densities had been increased by salt rejected by newly forming 
ice, was not able to mix vertically until density instabilities occurred in 
the surface region, thus inducing convective overturn. Since the presence 
of a strong surface halocline would require higher surface-salinity accumula- 
tions and hence permit more ice formation and associated salt rejection 
before this surface layer ivould become unstable, a strong halocline would 
be expected to inhibit convective overturn, both in frequency of occurrence 
and depth of penetration in the lead. Added to this insulating effect is 
the fact that, as surface salinities increase, surface freezing temperatures 
correspondingly decrease such that freezing rates are proportionately lowered. 

The role played by the vertical temperature gradient in lead thermodynamics 
was strongly dependent upon the existing vertical salinity gradient. Since 
density changes in sea water at temperatures near the freezing point are 
almost totally a function of salinity changes, the insulating qualities of 
the thermocline are undefinable when treated independently of the halocline. 
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However, in situations wherein weak vertical salinity gradients (poor in- 
sulation) and strong vertical temperature gradients coexisted, convective 
overturn resulted in relatively large quantities of heat being transported 
to the surface. Freezing and associated salt rejection were thus inhibited 
and surface melting was induced. In this manner, the temperature stratifi- 
cation was superimposed on the primary control exerted by the vertical 
salinity gradient. 

The effect of horizontal current magnitude on lead stability appeared, 
as in the case of the thermocline, to be secondary to the effect of the 
halocline. Strong currents will advect surface salt accumulations out of 
the lead relatively quickly, and will in general allow less time for lead 
thermodynamics to affect a given parcel of water advectively transiting 
the lead's width. This allows less change in the parcel's physical charac- 
teristics, Through these effects, the current also acts as an insulating 
agent. 

Ice growth rates were clearly dependent upon combinations of the above 
three variable parameters. Neglecting special ice profiles across the lead 
temporarily and examining only the maximum ice thicknesses occurring on the 
lead's surface as a function of time, the highest growth rates were exhibited 
for cases where a strong current, a weak thermocline, and a weak halocline 
occurred together. Conversely, lowest growth rates resulted from the co- 
existence of a strong current, a strong thermocline, and a strong halocline. 
Very little difference however, was noted in the high growth-rate case 
when currents were weak instead of strong, nor in the low growth-rate case 
when currents were weak or when currents were weak and the halocline was 
v;eak. Analysis of these results showed that high growth rates v;ere in all 
cases associated with v;eak vertical temperature gradients and low growth 
rates with strong vertical temperature gradients. Further, that current 



strength was relatively unimportant in changing growth rates when either 
strong haloclines were present, or when weak temperature gradients together 
with weak salinity gradients were present. 

The effect of the near-surface vertical temperature gradient appears 
then to dominate in determining surface freezing rates. Strong gradients 
resulted in lower rates of ice formation than in the case of weak gradients. 
Vertical transport of heat from depth in the lead by both eddy diffusion and 
convective overturn is then modified by both current strength and the near- 
surface halocline acting as insulating agents. 

Ice thickness profiles across the lead's width varied considerably with 
different combinations of the three basic variable parameters. Uniform or 
nearly-uniform ice thicknesses across the lead were in all cases associated 
with strong near-surface haloclines. Variations occurred where weak halo- 
clines existed, orowing larger when weak lialoclines and strong thermoclines 
occurred together. The greatest surface irregularity, both by amplitude of 
variation as well as by spacial frequency of variation, was found when weak 
currents, weak haloclines, and strong thermoclines occurred together. 

Analysis of these cases indicated that the halocline was the dominant 
variable in determining uniformity of ice thickness across the lead through 
insulation of the surface from the heat at depth in the lead. When convective 
overturn occurred as a result of increases in surface salinity, surface freezing 
rates v/ere altered at the point of convective overturn and produced surface 
irregularities. Since rejected salt is being continuously advected down- 
stream, there will exist horizontal salinity gradients which, when exceeding 
the magnitudes of vertical salinity gradients, will induce convective overturn 
and resultant surface melting at that point. Downstream of this point. 
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surface waters are diluted due to the upstream melting, thus experiencing 
increased freezing rates since their freezing point has effectively been 
raised. 

The near surface thermocline appeared to affect surface-ice uniformity 
through determining the amount of heat transported to the surface during 
convective overturn, and thus the amount of melting and dilution occurring. 
Finally, current strength appeared to act analogously to the halocline as an 
insulating agent, though less effectively. By flushing surface-salt accumu- 
lations out of the lead more rapidly, a strong current appears to inhibit 
convective overturn and hence variations in ice uniformity. By decreasing 
the amount of change a given parcel of water experiences over a given hori- 
zontal distance in the lead, a strong current appears to produce less 
intensive but more extensive lee effects from surface dynamics. 

Surface heat-loss densities varied across the lead's surface in response 

to ice thickness profiles. Since open-water sensible, latent, and radiative 

heat-flux densities used in this work totalled 2.00 x 10 ^ cal/cm^-sec as 

— ^ 

compared to a corresponding value over the 200 cm of ice of 1.15 x 10 
cal/cm^-sec, resultant heat-loss densities v/ere considerably larger over 
regions of the lead surface covered with relatively thinner ice. Uniformity 
of surface heat-loss densities over the lead's surface may then be directly 
related to the uniformity of the ice profile over the surface, while relative 
magnitudes of heat-loss densities may be inversely related to relative ice 
thicknesses . 

The extreme case of lead instability and non-uniformity (weak current, 
strong thermocline, and v;eak halocline) requires special comment. A high 
mean surface heat-loss density characterized this case, which was expected as 
a result of its highly non-uniform surface-ice profile. However, the lack 



i 







of evidence of these highly non-uniform ice covers occurring in nature in the 
absence of mechanical deformation forces prompted further analysis. By 
studying model output fields of runs 2, 6, 9, 10, and 14 (particularly the 
horizontal density gradients generated in the lead), an upstream pressure 
force was found to exist in the vicinity of the source boundary which was 
capable of performing on the order of ten ergs of vwrk against a unit mass 
of water entering the lead at the source boundary. Since a thermodynamically 
generated force of this type would completely alter the hydrodynamics 
occurring in the lead for low current magnitudes, the model appears to be non- 
representative of actual lead thermodynamics when using very low current 
magnitudes with strong thermoclines and weak haloclines. Stronger currents, 
on the other hand, had generally weaker horizontal density gradients associated 
with them and were also proportionately less affected by these pressure forces 
since their capacity to do work against the horizontal pressure force in- 
creases as the square of their magnitude. 

Results of run 14 proved to be of particular interest through the in- 
sights they afforded to the process of convective overturn. Since no current 
was present, there was no horizontal coupling between model grid points. As 
a result, although freezing progressed uniformly over the surface of the lead, 
strong differential horizontal density gradients developed along the 'source' 
boundary as cooling progressed. Of prime interest however, was the graphic 
illustration of the role played by the current as a convective-overturn 
scaling factor. With zero current, a strong temperature gradient and a 
weak salinity gradient, freezing- induced convective overturn penetrated to 
a depth of 35 meters roughly 18 hours after lead generation. Without the 
flushing action of a horizontal current or the insulating effect of a strong 
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halocline, the effects of surface thermodynamics were felt much deeper in 
the lead. It appears probable then, that in cases v^herein weak currents and 
weak haloclines exist in open leads (marginal vertical stability), surface 
freezing will induce, through convective overturn and vertical mixing, sub- 
stantial vertical and horizontal gradients in both temperature and salinity. 
As a result, significant fluxes of heat, salt, and momentum would be genera- 
ted in the waters of the lead, and sophistication beyond that possessed by 
the work at hand would be required to accurately represent the resultant 
lead hydrodynamics. 

Sunmarizing the analysis of results from integrating the open-lead model 
over the experimental ranges of variables, the near-surface vertical salinity 
gradient appears to exert the dominant force in maintaining vertical density 
stability in the water column through its insulating effect. By inhibiting 
convective overturn, a strong halocline was associated with horizontal ly- 
uniform ice thickness profiles and consequently uniform surface heat-flux 
densities. Current strength and the vertical temperature gradient appeared 

to be secondary as insulating agents. 

The near-surface temperature gradient appeared to be the dominant 
variable in controlling the rate of ice formation since its strength 
directly influences the amount of heat transported to the ice-water boun- 
dary from depth by both eddy diffusion and convective overturn. In this 
case, the halocline and the current strength appeared secondary m their 
effects on the freezing rate. 

Mean surface heat-loss densities were higher for more non-uniform 
surface ice thicknesses, which in turn may be inversely related to the 
strength of the near-surface halocline. This conclusion could prove to be 
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of significance to critical examinations of the heat budget over large 
areas of the arctic in which significant regional variations in the near- 
surface halocline might reasonably be expected. 

Representative freezing rates, ice-thickness profiles, and surface 
heat-loss densities found in integrating the model using real-case variables 
agreed well with values in the literature. Due to the similarity of current 
strengths, near-surface thermoclines , and near-surface haloclines, the fore- 
going relations v/ere generally not discernable, although the effect of down- 
stream salt accumulation was noticed in slightly decreasing dov^nstream ice 
thicknesses . 

Finally, it is concluded that, with the exception of cases in which 
weak horizontal currents occur together with strong thermoclines and vjeak 
haloclines, the model presented herein provides an accurate physical 
representation of the thermodynamics occurring in an open lead in central 
arctic sea ice over the experimental range of variables. 
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VI . RECOMMENDATIONS FOR FURTHER RESEARCH 



Although the surface boundary of the model could be considerably more 
sophisticated in its treatment of energy fluxes and transport processes, 
models already exist which have this sophistication, notably that of Maykut 
and Untersteiner [1969, 1971]. Coupling the model at hand to their model 
would therefore seem a more reasonable approach to refinement of the surface 
boundary in order to model air~ice“Sea interactions and mass~energy balances. 

The hydrodynamics occurring in a refreezing lead however, are still not 
accurately modelled, nor are they fully understood despite the insights 
afforded by the research performed. Extended development of this work 
would promote further conceptual understanding of these processes as well 
as provide a more physical ly-accurate model. Incorporation of accelerations 
resulting from thermodynamically-generated horizontal pressure gradients, 
as well as the coriolis acceleration, would considerably refine the model. 
Additional refinement would be gained by removing previously discussed 
limitations on: 1) the minimum vertical salinity gradient able to be 

accurately treated, and 2) the preclusion of downstream advection of warmer- 
then-freezing-point waters under an ice cover. 

In summary then, further research could profitably be applied to the 
work at hand to further refine the hydrodynamic accuracy of the model, as 
well as to couple it to larger-scale models of air-ice-sea interactions and 
mass-energy balance predictive schemes. 



APPENDIX A 



COMPUTER PROGRAM DESCRIPTION 

The computer program is written in FORTRAN IV and v/as used with the 
IBM 360/67 computer system at the W.R. Church Computer Center, Naval 
Postgraduate School. The program has been generalized such that sea 
water temperature and salinity profiles used as initial and left boundary 
fields are treated as input data. A design feature of the program allows 
the use of an internally variable time step to facilitate adjusting the 
model's temporal resolution. Additionally, the model allows variable grid 
spacing as well as variable lead dimensions, requiring changes only in the 
DIMENSION, COMMON and FORMAT statements. Although all currents used in 
the model were constant in time and space, the model allows use of both a 
depth-dependent and time-dependent horizontal current, again internally 
variable. With only minor modifications, the program could be easily ex- 
tended to include a bottom-boundary heat and salt flux, more sophisticated 
treatment of air-sea interactions at the surface boundary, and feedback of 
accelerations, resulting from horizontal density gradients generated by lead 
thermodynamics, to the motion of water through the lead. 

The program is divided into two basic sections: the main program and 

subroutine OVRTRN. The main program is further subdivided into seven parts: 
problem initialization, data input, computation of time-dependent temperature 
and salinity fields, ice-accretion/melting and salt rejection/dilution, with 
vertical stability check and computation of resultant heat flux densities, 
filling temperature and salinity arrays with resultant values, computation 
of applied surface heat flux densities and conversion to an integrated surface 



70 



temperature gradient, and writing output. Subroutine OVRIRfi accomplishes 
vertical convective overturn and mixing of unstable regions in the v/ater 
column. 

FORTRAN IV symbols for primary program parameters are defined below: 

AH Downstream-advected heat-flux density (cal /cm -sec) 

2 

AHF Applied surface heat-flux density (cal /cm -sec) 

2 

BTC Sea ice thermal conductivity constant (cal-cm^/g-sec) 

2 

CAH Cumulative downstream-advected heat density (cal /cm ) 

fy 

CAHF Cumulative applied surface heat density (cal/cm^) 

CHC Cumulative subsurface cooling heat density (cal /cm ) 

CHF Cumulative applied surface heat density less cumulative 
downstream-advected heat density (cal/cm^) 

2 

CHIF Cumulative heat density of surface-ice formation (cal/cm'^) 
Cl Thermal conductivity of sea ice (cal/cm.-sec-°C) 

CL Sea-water chlorinity (g/kg) 

CO Thermal conductivity of pure ice (cal/cm-sec-OC) 

CP Specific heat of sea water (cal/g-°C) 

DI Rate of surface ice formation (cm^/cm^-sec) 

DSALT Salt rejected per cm^ of ice formed (g/kg) 

DTAU Time step (sec) 

DTEMP Air-v/ater temperature difference (°C) 

DX Horizontal grid spacing (cm) 

DZ Vertical grid spacing (cm) 

q 2 

EPS Allowable rate-of-freezing error (cm-^/cm -sec) 

O O 

ERR Allowable ice thickness error (cmvcm"^) 

GRADI Vertical ice temperature gradient (oC/cm) 
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H 

HC 

HF 

HIF 

HO 

HS 

MHl 

MVl 

MV HI 

QI 

QL 

QR 

QS 

RHO 

RHOI 

S 

SALI 

SI 

SIGT 

SSCHG 



Model depth extent (cm) 

p 

Subsurfsce cooling heat-flux density (cal/cm -sec) 

Applied surface heat-flux density less downstream-advected 
heat flux density (cal /cm'^-sec) 

Heat-flux density of surface-ice formation (cal/cm^-sec) 
Heat of formation of pure ice (cal/g) 

Heat of formation of sea ice (cal/g) 

Number of grid columns 
Number of grid rows 
Total number of grid points 

Heat-flux density through sea ice to the atmosphere 
(cal/cm^-sec) 

Open-water latent heat-flux density to the atmosphere 
(cal/cm^-sec) 

Open-water radiative heat-flux density to the atmosphere 
(cal /cm^-sec) 

0pen-v7ater sensible heat-flux density to the atmosphere 

p 

(cal/cm -sec) 

p 

Density of sea water (g/cm ) 

Density of sea ice (g/cm^) 

Sea-water salinity (g/kg) 

Sea-ice salinity (g/kg) 

Sea-ice salinity (g/cm^) 

Sigma-t 

Salinity change in surface layer of depth DZ (g/kg) 



SSFLUX Salinity change in one cm^ of sea water at the ice-water 
boundary (g/kg) 

STGD Total applied surface temperature gradient at the water's 
surface (°C/cm) 

T Sea-water temperature (°C) 

Tau Time (sec) 

TF Sea-water freezing temperature (°C) 

TGRADL Surface temperature gradient due to latent heat-flux 
density (°C/cm) 

TGRADR Surface temperature gradient due to radiative heat-flux 
density (°C/cm) 

TGRADS Surface temperature gradient due to sensible heat- flux 
density (oC/cm) 

TI Temperature of lower ice surface (°C) 

U Horizontal current magnitude (cm/sec) 

VEC Vertical eddy conductivity (cal/cm-sec-°C) 

VSK Vertical salt eddy diffusivity (cm'^/sec) 

VTK Vertical temperature eddy diffusivity (cm /sec) 

W Lead width (cm) 

Z Depth (cm) 

ZI Surface ice thickness (cm^/cm^) 

A descriptive flow diagram of the program is shown in Figure 16. 
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Figure 16. Descriptive Flow Diagram of Program. 
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Figure 16 (continued). 
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ATI (B) = 1.-DTAU=? (U( B) / DX+2.=^VTK/0Z**2 ) 
ASKB ) = le-OTAU'--'(U( 8)/DX+2.=«'VSK/DZ*=4'2) 
AT2(B)=0TAU*U( B)/DX 
AS2(B)=0TAU*U( B)/OX 
AT3( D) =DvVTK/DZ--2 
AS3(D) =D*VSK/0Z'*“2 
C(A)={ A-0.03C) /lo8050 

C‘/\ / A \ A /. “7AQi!cA- 



DE(A) = loC’+{ A/10D3. ) 

F( A) =0«0-OoC966^A + 0o0000052=<'Av>i=3 ) 
R(C)=OSALT*C 
ERR(D) =D^=EPS 
BTC=0* 2 8 
C0=0,0C486 
DTAU=9C<,0 
DTEMP=25o 
DX=10Cvo0 
DZ=5C0«0 
EPS=1.0E-8 
H= 5000.0 
H0=79.77 
QL=3.0e-3 
QR=2.0E-3 
QS = 1. 5E-2 

RH0I=D.9i 
SALi=^8oO 

si=rhoi^sali/icD!;.o 

VSK=l0o0 
VTK=10,0 
W= 20000.0 
MH=IFIX( W/DX) 

HV=IFI X(H/DZ) 

MV1=MV+1 



r 

i o i. .132A) ) ) 
,000001324«A=^^=<‘3 



NF1=MK+1 

READ( 5, 100) IRUM 
IF( IRUN.EQ.99) GO TO 1000 
READ (5, 15";) (Tl'iEWd ) ,I = 1,MV1) 
REA0(5, 150) (SNEW( I ) ,I=ltMVl) 
CL=C( SNEW( 1 ) ) 

TF( 1)=F(CL) 

TNEW( 1 )=TF ( 1 ) 



N = iMVl 

DO 3 J=2,MH1 
DO 2 K=1,MV1 
I=M+K 

TNEW( I )=TNEW(K) 
T( 1 )=TMEH(K) 
SNEWd )=SIMEW(K) 
S( I )=SNEW(K) 

2 CONTINUE 
M=M+MV1 

3 CONTINUE 

DO 4 1=1, MVl 
CL=C( SNEW( I ) ) 
SIG0=SC (CL ) 
ATMP=AT (TNEW( I ) ) 
BTMP = B1 (TNEW( I ) ) 
SMT=ST(TN5W( I) ) 
SIGT( 1 ) = SGT(SVIT ) 
T( I )=TNEW( I ) 
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S( I )=SN:EW( I ) 

4 CONTINUE 
DO 5 

AH( I )=O.C 

AHF( I ) =OoO 
CAH( I )=0.0 
CAHF( I )=C,0 
CHC(I)=0.0 
CHF( I)=0.0 
CHIF( I )=0.O 
DI (I )=OoO 
HC( I )=0.0 
HF( I )=OoO 
HIF( I )=OoO 
ZI( I )=0»0 

5 CONTINUE 
J=MHl-( MH/A) 

DO 6 I=J,MH1 
ZI(I)=20Co0 

6 CONTINUE 
TAU=0.0 
GO TO 201 

101 TAU=TAU+OTAU 
DO 10 2 1=2, MHl 
AHF( I)=HUF( I ) 

CAHF( I)=CAHF( I ) + PUF( I) 

102 CONTINUE 
J=MV1+1 

DO 103 1=2, MH 

Z=0.0 

J1=J-MV1 

TNEWU ) =AT1 ( Z) ’^T ( J ) +AT2 ( Z ) *T (J 1 ) 

1 +AT3( DTAU>-2o*(T(J+1)~DZ*STGD(I)I 
SNEW(J)=AS1(Z)”S(J)+AS2(Z)>^S(J1) 

1 + AS3( OTAU) =!‘2o=!'S( J + 1) 

Z = H 

J=J+MV 
Jl= J-MVl 

TNEWU )=AT1(Z)»T(J) +AT 2 ( Z)^T ( J1 ) 

1 + AT3( DTAUI-Zo^i'TI J- 1 ) 

SNEW( J) = AS1(Z)*S(J )+AS2( Ziyi'S! Jl) 

1+AS3( DTAU)*2o*S( J-1 ) 

J=J + 1 

103 CONTINUE 

Z=OoO 
K=MVH1-HV 
J — K M V 1 

TNEw'( K) = AT1(Z)*T(K) -(-AT2 ( Z ) ’i^T ( J ) + AT3 ( OTAU ) ^2. 
1 (T (K + 1 )-DZ^''STGD( MHl ) ) 

SNEWIK )=AS1 ( Z) vS (K ) +AS2( Z)^S (J ) +A$3( OTAU ) -2. 
Z = H 

K=MVH1 

TNEW?kUaT1( Z)=«'T(K)+AT2(Z)=^T(J)+AT3<DTAU)+2o 
SNEW( K ) = ASl ( Z ) -S ( K ) + AS2 ( Z)^S ( J ) + AS3{ DTAU ) •i'Zo 
DO 105 J=2,MH 
Z=DZ 

Kl=( ( J-1 )*MV1) +2 
K2=( J^MVl ) -1 
DO 104 K=K1,K2 

TNEwTk'(=ATK Z) ’f'TIKI+ATZI Z)*T (J1 )+AT3(0TAUI* 

HNEmKUIsitzUsiKI+ASZIZI^SIJD+ASSIOTAU)- 
1(S(K-1)+S(K+1) ) 

Z=Z+DZ 

104 CONTINUE 

105 CONTINUE 
K1=1+(MVH1-MV) 

K2=MVH1-1 

Z = DZ 

DO ?.C6 K=K1,K2 



^;=S(K+1) 



•-^T(K-l ) 
4S(K-1) 
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J=K-M\/i 

TNEVJ(k1 = AT1(Z)*T(K)+AT2(Z)^T ( J) + AT3(DTAU)* 

^SNEwTk)=AS1(zUs(K)+AS 2(Z)=«=S (J)+AS3(0TAU 
1(S(K-1 )+S(K+l) ) 

Z=Z+DZ 

106 CONTINUE 
201 J=1 

DO 213 I=2,MH1 
J=J+MV1 
J3=J+MV 
MELT=G 

HIF(I)=OoC , 

THF=DZ/2o*'MT(J )-TNEVi( J ) ) 

2C2 CL=C(SNEW( J) ) 

TF( I )=F(CL) 

TTEMP=TNEW( J ) 

IFdTEWPoLEoTF ( I ) )TTEMP=TF( I ) 

IF(TTENp!gtItF ( ! )IaNd! (H ( I ) oLE. ER5^{ DTAU) ) )TIF=C*0 
CP=SH(SNEW( J) ) 

SIGC=SC(CL ) 

ATMP=AT (TTEMP) 

BTMP=BT( TTEMP) 

SMT=ST(TTEWP) 

SIGT( J)=SGT(SMT ) 

RH0=DE{ SIGT( J ) ) 

ZKW=P.HO''CP/DTAU 
SRAT10=SAL I/SMEW( J ) 

HS=HC*( 1.0-SRATI 0) 

HICE=(TIFvZKW) /( PH01*HS) 

GO TO 20^ 

203 HICE=-( ZI( I )/DTAU) 

TIF=HlCE*RhCI^HS/ZKW 

204 G1=-(DTAU^=H1CE ) 

IF(GUGT.ZI( I) ,'GC TO 203 
0SALT=KH0d (SNEWv J )-SALl ) 

ZICE^Dl AU*HiCE 
SSFLUX=R(ZICE) 

SSCHG=SSFLUX/DZ 
IF(NELT«EQ«r)GO TO 206 
SDIST=SSCHG/JS A 

Kl=J+( JSA-1 ) 

DO 205 K=J,K1 

SNEW( K)=SNEW(K)+SDIST 

205 CONTINUE 
GO TO 207 

206 SNEW( J)=SNEW( J )+SSCHG 

207 TNEW( J)=TMEW( J)+(TIF*2e/DZ) 

iM in I ilLTtERRVoTAufoANO* ( H IC E» L E «0 «0 ) )ZI ( I )=QoO 
HIF( I)=HIF( I )+TIF-ZKW 
IF(ABS(HICEI,GTeEPS)GO TO 2o2 
CL=C(SNEW( J) ) 

TF( I)=F(CL) 

DO 208 K=J,J3 
CL=C(SNEW(K) ) 

SIGC=S0{CL) 

ATMP=AT(TNEW(K) ) 

BTMP=BT (TNEW(K) ) 

SMT=ST( TNEW( K) ) 

SIGT(K)=SGT(SMT) 

208 CONTINUE 
IF(MELT.EQo 1 )G0 TO 212 
J2=J3+1 

209 DO 210 K=1,MV , 

IF(SIGT( J2-K).GE«SIGT( J3-K) )G0 TO 210 
L=J2-K 

CALL OVRTRN (J,L) 

GO TO 209 

210 CONTINUE 
CL=C(SNEW( J) ) 
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TF( I)=F(CL) 

IF(TNEW( J) ,LE.TF( I ) )G0 TO 212 
IF( Zi ( I ) .LtoERR( DTAU) ) GO TO 212 
MELT=1 

DO 211 K=1 ,MV 

IF(SNEW( J+K)«EQoSNEW(J) )GO TO 211 



JSA=K 
GO TO 202 

211 CONTINUE 
JSA=MV1 
GO TO 202 

212 OK I)=HIF( I )/( RHOI*HS) 

CHIF( I )=CHIF( I)+HIF( D^S'DTAU 
HF( I)=THF>!‘ZKW 

CHF( I) =CHF( I) + HF( I )#DTAU 
AH( I ) = AHF( I )-HF( I ) 

CAH( I)=CAH( I )+DTAU-AH( I ) 
HC(I)=HF( I)-HIF( I ) 

CHC( I )=CHF( I )-CHIF ( I ) 

213 CONTINUE 
I 1=MV1+1 

DO 301 1=11, NVHl 
T( I )=TNEW( I ) 

S( I ) = SNEW( I ) 

301 CONTINUE 
J=MV1+1 

DO 403 I=2,MH1 
RHO=OE(SIGT( J) ) 

CP=SH(S( J) ) 

VEO = VTK’'"RHQ^CP 

IF(ZI( I ).LEoERR(DTAU) )G0 TO 401 
TI=TF( I ) 

CI = C0+3TC=^SI/TI 
GRADI = DTEMP/ZI ( I ) 

QI=CI=<=GRA0I 
TGRApS = QT/'v'rc 
T GRADL=0 c'.i 
TGRADR=0o55^QR/VEC 



TGOW=QS/VEC 

IF(TGRAOSoGToTGOW)GO TO 401 
GO TO 402 

401 TGRADS=QS/VEC 

TGRADL = QL#EXP(-0o5--i'ZI ( I ) )/VEC 

TGRADR = (0. 35''<QR+0,65-^QR^^EXP(-C.5=i'ZI( I) ) )/VEC 

402 STGD( I ) = TGRADS+TGRADL+TGRADR 
HUF( I )=STGD( I I'^'VEC 

PUF( I ) = HUF( I I'^DTAU 
J=J+MV1 

403 CONTINUE 

IF(AM0D(TAU,3600o ). LT* l.E-4) GO TO 501 
GO TO 101 



501 Jl=l 

TIME=TAU/3600o 



J2=MVH1-MV 

WRITE (6,20D)TIME, IRUN 
WRITE(b,25C) (TF(K) ,K=1,MH1) 

DO 502 1=1, MVl 

WRITE! 6, BOOM T( J) , J = J1 , J2 , MVl ) 

WRITE (6, 35. ■ ) (S ( J > , J=J1, J2,MV1) 

WRITE (6, 400) (SIGT( J ) , J=J1 , J2,MV1 ) 

J1=J1+1 

J2=J2+1 

502 CONTINUE 
WRITE! 6,450 
DO 503 I=2,HH1 

WRITE! 6, 500) HI F ! I ) , CHI F !I ) , DI ! I ) , Z I ! I ) . HC ! I ) 
1HF(I ) ,CHF( I ) , AH! I ) , CAH ! I ) , AHF! I ) ,CAHF! I ) 

503 CONTINUE 
IF!TIME,GEo48. )G0 TO 1 
GO TO 101 

100 FORMAT! 12) 

150 FORMAT! 11F7» 2) 



CHC! I ) , 
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200 FORMAT (• 1* I 'PROF- ILE ( X-Z PLANE) GF LEAD AT TI«E = »,F5*1 
1 tlOX, 'RUN NUMBER* ,13) 

250 FORHAT( 'O' , • TF ' , 3 X , 21 ( 1 X , F 5* ?. ) ) 

300 FORMAT (•C*,'T',4X,21(1X,F5.2)) 

3 50 FORMAT (• • , ' S ' , 4X, 2 1 ( 1 X , F5 . 2 ) )- 
400 FQRMAT(* * , * ST • , 3X , 21 ( 1 X, F 5. 2 ) ) 

450 FORMAT(//, 'C , 3X , * H I F • , 7 X, • CHI F • , 6X, » DI ' , 8X , ' Z 1 * , 8X , 
I'HC* ,8X, 'CHC, 7X, *HF* ,8X,*CHF' ,7Xt*AH* taX, •CAH*,7X, 
2'AHF* ,7X, *CAHF • , /) 

500 FORMATC • , 1 PE9, 2 , 1 1 ( 1 X , E9o 2 ) ) 
lOCO STOP 
END 



SUBROUTINE CVRTRN (I,K) 

CCMMON TNEWI232) ,SNEU(232) ,SI6T(232) 

C(A) = ( A-C.<.'30)/lo8( :5M 

SQ( A)=-C,069+lc4708=i=A-0o001570*A**2+0o00C0398*A^*3 

AT(A) = A* (4o 786 7-^ 0 098 185*A+0c 00 1 0843*A=!=’i= 2 ) . 001 

BT ( A ) =A4 ( 1 8 o03 C-Oo 8 164=«'A+C' o 0 1667*A=i==<'2 ) ’i'Cc OCCOOl 

ST(A)=-( A-3o98 )**2/503« 570^( A+283* ) / ( A-i-67e 26 ) 

SGT(A)=A+(SIG0+0ol324)«(l,0-ATMP+(BTMP*( STGO-0.1324) ) ) 

K1=K-1 

J$=l 

L=0 

DO 1 M=I ,K1 
L = L + 1 

IF(SIGT(K-L)oGE.SIGT(K))GO TO 1 
Ll=l+( K-L) 

J$=2 
GO TO 2 

1 CONTINUE 

2 CTEHP=0.n 
CSALT=C,C 

IF( JS. EQ.l )L1=I 
DO 3 N=I.1,K 

Cl EMP=C I BMP'* i NEvi ( N ; 

CSALT=CSALT + SNEV,'( N) 

3 CONTINUE 
NN=1+(K-L1 ) 

ATEMP=CTEMP/FLOAT( NN) 

ASALT=CSALT/FLOAT(NN> 

CL=C(ASALT) 

S1G0=SC (CL ) 

ATMP=AT(ATEMP) 

BTMP=BT (ATEMP) 

SMT=ST( ATEMP) 

ASIGT=SGT( SMT) 

DO 4 M=L1,K 
TNEW(M)=ATEMP 
SNEW(M )=ASALT 
SIGT(M)=ASIGT 

4 CONTINUE 
RETURN 
END 
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